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Ferrogels, i.e. swollen polymer networks into which magnetic particles are immersed, can be 
considered as “smart materials” since their shape and elasticity can be controlled by an external 
magnetic field. Using molecular dynamics simulations on the coarse-grained level we study a ferrogel 
in which the magnetic particles act as the cross-linkers of the polymer network. In a homogeneous 
external magnetic field the direct coupling between the orientation of the magnetic moments and 
the polymers by means of covalent bonds gives rise to a deformation of the gel, independent of the 
interparticle dipole-dipole interaction. In this paper we report quantitative measurements of this 
deformation, and in particular, we investigate the gel’s elastic moduli and its magnetic response for 
two different connectivities of the network nodes. Our results demonstrate that these properties 
depend significantly on the topology of the polymer network. 


I. INTRODUCTION 

Polymer gels combined with magnetic nanoparticles 
are composite soft materials whose elastic properties can 
be controlled by external magnetic fields [? ]. In the 
last decade an active research in this area put forward 
several terms to characterise these materials: magnetoe¬ 
lastomers, ferrogels magnetic gels, etc. These materials 
gained attention both, in fundamental research and ap¬ 
plications, as they combine the intrinsic properties of a 
hydrogel with a relatively high sensitivity to magnetic 
fields inherent to magnetic fluids [? ]. Of particular im¬ 
portance is the potential to change the shape and me¬ 
chanical characteristics of a ferrogel by applying a rela¬ 
tively weak magnetic field [??????]. 

The ability to control the shape of a magnetic gels gave 
rise to a number of new applications. Among them are 
drug delivery [? ? ], cancer treatment strategies (hy¬ 
perthermia) [? ], actuation [? ? ] and transport [? ? 
]• 

On a microscopical level, three ingredients need to be 
considered: the magnetic fluid, the polymer network, and 
the coupling between them. First, magnetic fluids (fer- 
rofluids) have been available for decades. They have been 
studied extensively experimentally and theoretically, and 
are used in technical and medical applications[? ? ? ]. 
The magnetic nanoparticles typically have a size on the 
order of 10 nm and hence consist of a single ferromag¬ 
netic domain. This implies that they carry a permanent 
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magnetic moment. Due to Brownian motion, the mag¬ 
netic particles in a ferrofluid are randomly oriented and 
the fluid as a whole does not exhibit a net magnetic mo¬ 
ment. This changes, once an external field is applied, 
and the magnetic moments in the fluid can align to it. 
There are two mechanisms by which this can occur [? ]: 
the magnetic moments of the particles can re-align in¬ 
ternally, without a rotation of the particle as a whole 
(Neel mechanism), or the particles as a whole can rotate 
to align their moments (Brownian mechanism). In the 
case of magnetic gels, it is of particular importance to 
know which mechanism is present, because only for the 
Brownian mechanism a direct coupling between the ori¬ 
entation of the magnetic moments and the polymers is 
possible. Moreover one can also study ferrogels, where 
the magnetic nanoparticles have a finite uniaxial mag¬ 
netic anisotropy [? ]. 

Looking at the second ingredient, hydrogels, i.e. water- 
filled polymer networks, are common in everyday appli¬ 
cations such as contact-lenses, food and medical items [? 
? ? ? ]. They can be synthesized from a wide range 
of polymers and their chemical and physical properties 
can be tailored to suite the applications. For instance, 
by varying the degree of cross-linking between the poly¬ 
mers, their stiffness can be controlled. Also, hydrogels 
can be made sensitive to environmental factors such as 
temperature and ph-value. In the context of magnetic 
gels, it is particularly important to distinguish between 
the type of cross-linking between the polymers. This can 
be either due to comparatively weak interactions such as 
van der Wals forces and hydrogen bonds, or it can be due 
to covalent bonds. While the links between the polymers 
are dynamically changing in the first case, in the latter 
case, the bonds are stable and the network topology does 
not change over time. Simulations of hydrogels have been 
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performed, but the literature as by far not as extensive 
as for polymer solutions [? ] 

The third thing to be considered when studying a fer- 
rogel is the coupling between the polymers and the mag¬ 
netic particles. One possibility is a loose coupling by 
means of hydrogen bonding or van der Wals interactions. 
In this case, the magnetic particles can exert a stress 
on the polymer matrix as they move, but the rotational 
degree of freedom of the nanoparticles is not directly cou¬ 
pled to the polymers. A rotational coupling between the 
nanoparticles and the polymers can be achieved by bind¬ 
ing the polymers covalently to specific spots on the sur¬ 
face of the magnetic nanoparticles[? ? ? ]. When the 
particles rotate due to an external field, the polymers 
wrap around the particles. This, in turn creates a stress 
on the polymer chain which leads to a shrinking of the 
gel. A sketch illustrating this mechanism can be found in 
Fig-lU In this case, the magnetic particles act as actual 
cross-linkers of the network[? ]. 

In this paper, we consider such a magnetic nanoparticle 
cross-linked gel, which due to the covalent bonds has a 
network topology which is constant over time. The mag¬ 
netic moments are assumed to relax via the Brownian 
mechanism, i.e., by a rotation of the particle as a whole. 

Since many applications of ferrogels are connected to 
actuation in an external magnetic field, it is crucial to un¬ 
derstand the microscopic mechanisms by which the gels 
are deformed. First, in a field gradient, the magnetic par¬ 
ticles tend to accumulate in regions with a stronger field, 
due to magnetostriction [? ]. Second, in a homogeneous 
field, a deformation can occur due to a change of the 
interactions between the magnetic particles as they are 
aligning in an external field. The dipole moments in the 
gel system are randomly oriented in the field-free case, 
and there is on average only a weak interaction between 
them, which is not strong enough to overcome the elas¬ 
ticity of the polymer matrix. The interaction between 
magnetic particles can become significant, however, as 
their moments are aligned in a field. Due to this change 
in interaction, the particles rearrange and thereby deform 
the matrix. This is the microscopic equivalent of the fact 
that a homogeneously magnetized elastic medium elon¬ 
gates in field direction in order to reduce its demagne¬ 
tization field[? ? ? ]. This mechanism was observed, 
e.g., in Ref. [? ]. It is noteworthy that the type of local 
particle configuration influences the response to a fleld[? 
? ]. Finally, a third deformation mechanism arises from 
a direct coupling of the orientation of the magnetic mo¬ 
ments to the polymers as described above. A detailed 
study on the interaction between two magnetic particles 
connected by a polymer chain in this way can be found 
in [? ]. 

In our earlier work[? ], we presented two gel models 
in two dimensions, which deform by the second and third 



FIG. 1. Sketch of magnetic nanoparticles in a gel. Upper 
row (a) and b)): 2D model, left without an applied external 
field, right with an applied external field, as marked with the 
red arrow. One can see four polymer chains rigidly bound 
to specific spots (light blue, virtual sites) on the surface of 
a magnetic nanoparticle. When a field is applied and the 
magnetic particle aligns its dipole to it, the polymer chains 
wrap around the particle. This exerts a stress on the chain, 
which leads to a shrinking of the ferrogel (isotropic in 2D). 
Lower row (c) and d)) 3D model, left without an external 
field, right with an external field applied. In sketch c) a 4-fold 
(diamond) and 6-fold (cubic) crosslinker is shown. In sketch 
d), an external magnetic field H is applied within the plane, as 
shown with the red arrow. Notice that only chains within the 
plane are subjected to the deformation, the polymer chains 
perpendicular to the plane are not affected. This leads to an 
anisotropic shrinking of the gel in 3D. 


mechanism, respectively. For the case of a gel cross-linked 
by magnetic particles (model II in Ref. [? ]), an isotropic 
shrinking can be observed. This is, because in 2D, there 
is only one rotation axis, namely the one perpendicular to 
the model plane. Once the magnetic particle is rotated by 
the held, all chains attached to it receive the same amount 
of stress. In Ref. [? ], however, we have shown that the 
shrinking is no longer isotropic in three dimensions. This 
can be explained by the effect that the rotation always 
takes place around an axis which is perpendicular to both, 
the magnetic moment prior to alignment and the direc¬ 
tion of the external magnetic held. Polymer chains that 
are attached within a plane spanned by these two vec¬ 
tors and that are going through the particle’s center fol- 



3 


low the full rotation of the magnetic particle. Chains 
attached parallel to the rotation axis, on the other hand, 
are not affected by the particle’s rotation. Prior to the 
alignment by a field, the magnetic moments of different 
particles are random, but the rotation axis around which 
the alignment takes place is always perpendicular to the 
external field. Hence, there is less deformation in the di¬ 
rections perpendicular to the field. A sketch explaining 
this mechanism can be found in Fig.[2 

In contrast to our previous works, we present here a 
detailed quantitative study of the deformation, elasticity 
and magnetic response of a three-dimensional gel cross- 
linked by magnetic particles, which has been investigated 
in depth in the PhD work of the first author [? ]. We con¬ 
sider two topologies of networks, a simple cubic (SC) net¬ 
work with six-fold connectivity and a diamond topology 
with fourfold connectivity (DC). The paper is structured 
as follows: In Sec.[ll| we describe the simulation technique 
used, and then introduce the model. Thereafter, the equi¬ 
librium degree of swelling of the gel is obtained (Sec. in. 
In Sec.[lV| the gel’s deformation in a field is measured. 
Sec. [V| deals with the measurement of the system’s elastic 
constants and in Sec. VT the magnetic response is exam¬ 
ined. We conclude the paper with a summary. 


II. SIMULATION TECHNIQUE 

In this paper, we study the deformation of three- 
dimensional magnetic gels cross-linked by magnetic 
nanoparticles in a magnetic field. This is done by means 
of molecular dynamics simulations on a coarse-grained 
level, using the ESPResSo software[? ? ]: rather than 
considering individual atoms as degrees of freedom, the 
smallest units considered in the simulations are the units 
of polymer that roughly match one Kuhn length , and 
the magnetic nanoparticles. In this way, the number of 
degrees of freedom in the system is reduced to a point, 
at which it is possible to simulate several meshes of the 
gel in each Cartesian direction. The simulations are per¬ 
formed in the canonical NVT ensemble using a Langevin 
thermostat. For each translational degree of freedom, the 
equation of motion is 


TflX — yX F ^random: ( 1 ) 

where m denotes the mass, 7 the friction coefficient, F 
the force due to other particles, and Frandom a random 
force mimicking collisions with the surrounding solvent. 
According to the fluctuation-dissipation theorem, in or¬ 
der to maintain a thermal energy of ksT^ each random 
force component has to have a mean of zero and a vari¬ 
ance of 


( 2 ) 


Rotational degrees of freedom are treated in a similar 
way. In Eqn.[^ position, mass and force are replaced by 
orientation, inertial tensor and torque, respectively[? ]. 

The polymer chains are represented as bead-spring 
models, in which the beads are bound by harmonic po¬ 
tentials 

Fharmonic('^) ~ 

where k is the spring constant and vq the equilibrium 
distance. 

Both, the beads forming the polymer chains and the 
magnetic particles, are represented as soft spheres inter¬ 
acting via a purely repulsive Lennard-Jones interaction, 
the so-called WCA-potential[? ]. Its functional form is 
given by 

UwcA = 

[0 r >rc 

where e controls the energy scale of the potential, a is the 
sum of the radii of the interacting beads, and Vc = 2^/^cr 
denotes the cut-off distance after which the potential is 
zero. 

The anchoring point of the covalent bonds between 
the magnetic nanoparticle and the polymer chains at¬ 
tached to it is modelled using the virtual site feature of 
ESPResSo [? ]:. Virtual sites are particles whose position 
is not determined by integrating an equation of motion, 
rather it is calculated from the position and orientation 
of other particles. In this particular case, a virtual site 
is placed immediately under the surface of the magnetic 
nanoparticle to form a binding site for the attached poly¬ 
mer. The latter is rigidly connected to the magnetic par¬ 
ticle and will follow both, its translational and rotational 
motion. The technical details can be found in the ap¬ 
pendix of Ref. [? ] as well as in [? ]. A sketch of a 
magnetic particle with the attached polymer chain can 
be found in Fig.[2 

The magnetic particles, forming the nodes of the net¬ 
work, are initially placed on a regular three-dimensional 
lattice. Then, the chains are attached to a specific point 
on the nodes’ surface. Therefore, when the node particles 
rotate, the chain ends have to follow. Before the chains 
are connected, the dipole moments of the node particles 
are randomized, to make the system isotropic. The sys¬ 
tem is studied taking into account only the dipole-field 
interaction, but not the dipole-dipole interaction. We do 
this for three reasons: first, the density of magnetic par¬ 
ticles is relatively low, so that the influence of the dipole- 
dipole interaction is rather weak. This was checked in 
the 2D case discussed in Ref. [? ]. Second, the inclusion 
of the dipole-dipole interaction could result in two defor¬ 
mation mechanisms occurring at the same time, which 
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FIG. 2. Snapshot of a gel sample at equilibrium swelling in the diamond cubic geometry with four chains attached to a 
node (left) and the simple cubic geometry with six chains attached to a node (right). Consisting of 60 beads, the chains are 
connected to specific spots on the surface of the magnetic particles. In the field-free case, depicted here, the magnetic moments 
are randomly aligned. 


would make it more difficult to attribute findings to a 
specific mechanism. Third, for particles placed on a reg¬ 
ular lattice, the kind of deformation due to dipole-dipole 
interactions depends on the lattice structure [? ]. In or¬ 
der to avoid artefacts, these interactions should therefore 
only be included in a more randomized system. 

In this work, two simple choices for an initial lattice 
are considered. These are the diamond cubic (DC) and 
simple cubic (SC) geometries. In the former case, four 
chain ends are connected to each node, while there are 
six in the latter case. By switching from the diamond 
cubic to the simple cubic structure, one can thus vary 
the degree of connectivity (and in consequence the elastic 
and magnetic response) of the network. 

In contrast to the 2D case, which we studied in model 
II of Ref. [? ], periodic boundary conditions are applied, 
in order to reduce surface effects, and to study a truly 
macroscopic gel. Snapshots of samples after equilibra¬ 
tion are shown in Fig.for the DC and SC lattice. The 
simulation parameters are as follows: the gel consists of 
Nn = 64 node particles as well as = 2Nn chains in 
the diamond cubic structure, or chains in the 

simple cubic structure. The chains consist of 60 or 80 
beads, with a diameter of dc = 1, whereas the node par¬ 


ticles have a diameter of = 10. The scale of the soft 
sphere interaction (Eqn.[^ is e = 10. The bonds between 
neighboring particles in a chain as well as between the 
ends of the chains and the virtual sites interact via a har¬ 
monic potential (Eqn.|^ with an equilibrium elongation 
of 26 ( 7 , coinciding with the cut-off of the WCA-potential. 
The spring constant of the harmonic potential is /c = 10. 

As the shape of the gel for the given parameters is not 
known a priori, an iterative procedure is applied to adjust 
the box size. We begin with an estimated shape. Then, 
iteratively, the stress is measured and the box is shrunk or 
expanded in each coordinate, according to the diagonal 
elements of the observed stress tensor. The orthogonal 
basis of the stress tensor coincides with the orientation 
of the simulation box. After the box shape is adjusted, a 
new simulation is performed. This is repeated, until the 
measured stress approaches zero with the desired toler¬ 
ance. When the gel can be assumed to be isotropic, i.e., 
when no external field is applied, an optimized procedure 
can be used, which will be explained in the next section. 
Even in the field free case, there will be a small anisotropy 
in an individual sample of the gel, because the sum of the 
randomly drawn initial dipole moments is never exactly 
zero. The stress is measured and averaged for several in- 
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stances of the system at any given set of parameters to 
compensate for this. In other words, the simulation is 
divided into the following steps: 

1. The magnetic node particles are placed on the lat¬ 
tice with random dipole orientations. The distance 
between the surfaces of adjacent nodes is such that 
the chains can be inserted in an elongated configu¬ 
ration. Hence, the resulting initial lattice constant 
is 

^ T (5) 

where and (Jc, are the Lennard Jones diameters 
of the node and chain particles, and Ic is the number 
of particles in a chain. 

2. Chains are connected to the surface of two adjacent 
nodes. Four chain ends are connected to a node in 
the diamond cubic case and six in the simple cubic 
case. 

3. After the network is cross-linked, it is scaled to the 
desired shape. This is done by adjusting the box 
size in 400 steps from the initial shape to the target 
shape. After each deformation step, the system is 
equilibrated with the Langevin thermostat for 100 
iterations of time step dtreshape = 0.0006 to disperse 
the energy introduced by the deformation. 

4. When the final shape is reached, the system is equi¬ 
librated further for 200 000 steps of size dt = 0.01, 
again using the Langevin thermostat. 

5. Finally, observables can be measured. The time 
step is again dt = 0.01 and stress measurements are 
taken every 20 time steps, whereas the magnetiza¬ 
tion is measured every 100 000 steps. In a typical 
simulation, more than 50 000 stress measurements 
are taken. This enormous amount of statistics is 
needed due to the large fluctuations of the mea¬ 
sured stress due to the softness of the material. 


III. EQUILIBRIUM SWELLING 

When a gel with specific topology, chain length, and 
interaction parameters is to be studied, the equilibrium 
swelling, and hence the box size, is not known a priori. 
Therefore, the first step in the analysis of the system is 
to determine this property in the field-free case. In prin¬ 
ciple, it is possible to do this by simulating in the NPT- 
ensemble using a barostat. However, this was found not 
to be efficient[? ]. Instead, a set of simulations in the 
A/'UT-ensemble is performed at different volumes. By 



FIG. 3. Stress-strain relations for both, the diamond cubic 
(filled symbols) and the simple cubic (open symbols) geom¬ 
etry. The x-axis is rescaled such that unity represents the 
equilibrium volume for the given geometry and chain length. 
From the respective slopes of the stress-strain relations at the 
point of equilibrium swelling, the elasticity of the samples can 
be determined. It can be seen that for both chain lengths 
shown, the samples in the diamond cubic geometry are softer 
than those in the simple cubic geometry. Additionally, longer 
chains make the gel more deformable. 

measuring the stress tensor and averaging over its di¬ 
agonal elements, the stress-strain relation for the gel is 
obtained: 

P{e) = \^Su, (6) 

where P(e) denotes the isotropic pressure at a certain 
strain e, and Sij denote the elements of the stress ten¬ 
sor. From this curve, the equilibrium strain, at which 
the stress is zero, can be estimated by linear interpola¬ 
tion. Taking the box length of the initial system to be 
unity, with completely stretched polymers, the following 
equilibrium swelling was obtained: for the diamond cubic 
lattice, values of 0.2790 and 0.2546 are found for chains 
with 60 and 80 beads, respectively. The values for the 
simple cubic structure are 0.3153 and 0.2742. It is con¬ 
venient, to express the size of the gels in terms of this 
equilibrium swelling. This will be done throughout the 
rest of the paper. 

In Fig.|^ the stress-strain curves for both, DC and SC 
geometries, are shown. From the slope of the curves at 
the equilibrium strain, the stiffness of the gel can be de¬ 
duced. The higher the slope, the stiffer the material is. 
It can be seen that the gel in the SC geometry is more 
stiff than in the DC geometry. Also, a decrease of the 
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FIG. 4. Stress measured parallel and perpendicular to the 
external field versus the strength of the external held. Data 
is shown for the diamond cubic (DC, full symbols) and simple 
cubic (SC, open symbols) geometry and for chain lengths (cl) 
of 60 and 80 beads per chain. It can be seen that in both 
geometries, the absolute value of the stress parallel to the 
held is larger than the value measured in the perpendicular 
direction. A much larger negative stress is observed in the 
simple cubic geometry than in the diamond cubic one. 


chain length results in a stiffer system. Both these obser¬ 
vations can be understood by noting the higher particle 
density for shorter chain lengths and for a larger number 
of chains. A quantitative study of the elastic constants 
will follow in Sec.lVl 


IV. ANISOTROPIC DEFORMATION 


helds are expressed in terms of the Langevin parameter 


a = 


fiomH 
ksT ’ 


(7) 


where m denotes the particles’ moment, /io the vacuum 
permittivity and H the external held. This parameter 
measures the ratio between the Zeeman energy and the 
thermal energy ksT. The measured stresses for all ap¬ 
plied external magnetic helds larger than zero are nega¬ 
tive, indicating that the gels would contract if the con¬ 
straint of the constant volume would be removed. The 
absolute value of the stress in the direction parallel to the 
held is signihcantly larger than the one measured in the 
perpendicular direction. It is also worth noting that the 
stress is signihcantly stronger in the SC geometry, where 
six chains are connected to a node as opposed to four in 
the DC geometry. This implies that two contradicting 
forces are at play. On the one hand, in the simple cu¬ 
bic structure, there are more chains attached to a node. 
This leads to a stronger wrapping effect and a stronger 
tendency of the system to deform in the held. On the 
other hand, from the slope at the equilibrium swelling in 
Fig.i it can be observed that the gel based on the simple 
cubic structure is stiffer than the one based on the dia¬ 
mond structure. From measurements of the deformation 
of the gel, it turns out that for this particular model, the 
elasticity dominates the deformation behaviour. 

To hnd the equilibrium shape of the gel for a given ex¬ 
ternal magnetic held, simulations in the WT-ensemble 
are executed iteratively. In each simulation, the stress 
is measured in the directions parallel and perpendicular 
to the external held. Then, the box size for the next 
simulation is altered according to these stress measure¬ 
ments. This procedure is repeated until the remaining 
stress approaches zero within the chosen error tolerance. 
The details of the procedure are as follows: 


In this section, the deformation of the gel under the 
inhuence of an external held will be discussed. In partic¬ 
ular, we will have a closer look at the anisotropic deforma¬ 
tion described in Ref. [? ]. As the hrst step, we measure 
the tendency of the gel to deform in the directions paral¬ 
lel and perpendicular to the held, due to the wrapping of 
the polymer chains around the magnetic particles (Fig.[^. 
In order to separate this wrapping effect from the elas¬ 
tic properties, we hrst measure the direction-dependent 
stress due to an external held in the non-deformed gel. 
Hence, the gel is simulated in the WT-ensemble at the 
equilibrium volume for the held free case as it was de¬ 
termined in the previous section. Then, a held is applied 
and the stress in the directions parallel and perpendicular 
to the held is measured. The results for chain lengths of 
60 and 80 beads, and for both, the DC and SC geometry, 
are shown in Fig.|^ Throughout this paper, magnetic 


• A shape of a simulation consists of the lengths of 
the simulation box parallel and perpendicular to the 
external held. 

• For each shape simulated, the stress measurement 
is averaged over 24 simulations to reduce the mea¬ 
surement error. 

• Initially, two sets of simulations are performed for 
a given external magnetic held. The gel samples in 
them have manually selected shapes, of which one 
is above and the other is below the expected equi¬ 
librium shape for the given held. For these shapes, 
the corresponding stresses in the direction parallel 
and perpendicular to the held are measured. 

• In all further iterations, a single new shape consist¬ 
ing of a box length parallel and perpendicular to 
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FIG. 5. Left: relative loss of volume of a gel versus the strength of the external field. Right: length of the gel in the directions 
parallel and perpendicular to an externally applied magnetic field versus the strength of that field. Full symbols denote the 
diamond cubic (DC) geometry, whereas open symbols are used for the simple cubic (SC) one. The chain length (cl) was 60 
(squares) and 80 (circles), respectively. It can be seen that the overall shrinkage is comparable in magnitude for all systems 
considered within the simulation error, which, however, is large. While all samples shrink significantly more in the direction 
parallel to the field, it is noteworthy that the gels in the simple cubic geometry shrink less in the direction parallel to the field, 
but more in the perpendicular direction than those in the diamond cubic geometry. This results in a similar change of volume 
in all systems, although the “shape “ of the deformation depends strongly on the network geometry. 


the field is generated. For the respective directions 
parallel and perpendicular to the field, a new box 
length is generated from the two ones examined in 
previous iterations which yielded the lowest corre¬ 
sponding absolute stress. From these box lengths 
and corresponding stresses, a linear stress-strain re¬ 
lation is constructed, which is used to extrapolate to 
the point where the stress would be zero according 
to this relation. Due to large statistical errors in the 
stress measurements, it is necessary to restrict the 
change in box length to approximately two percent 
per iteration, to keep the procedure stable. For this 
new shape, again, the stress is measured in the di¬ 
rections parallel and perpendicular to the external 
magnetic field. 

• The iteration is terminated, if the stress in both, the 
direction parallel and perpendicular to the field, is 
below 10“^. 

The iterative procedure provides satisfactory results but 
requires a large amount of computational resources and 
occasional manual interventions. In a future study, 
more sophisticated procedures, like genetic optimizations 
might be considered. As an alternative to the iterative 
procedure, it is also possible to simulate a set of shapes 
around the expected equilibrium shape, obtain the corre¬ 
sponding stresses, and fit a linear stress-strain model to 


the data. This is done in conjunction with the measure¬ 
ment of elastic constants in Sec. El 

To study the deformation of the gel in an external mag¬ 
netic field, the equilibrium shapes were determined for 
fields up to (a = 60 in terms of the Langevin parame¬ 
ter. Results for chain lengths of 60 and 80 beads are 
shown on the right hand side of Fig.for gels based on 
both, the diamond cubic and simple cubic network struc¬ 
ture. It can be seen that, as expected from the anisotropic 
stress measurements at constant volume (Fig. Q and as 
explained above, the gel deforms significantly more in the 
direction parallel to the external field than in the perpen¬ 
dicular direction. The deformation decreases when the 
chain length is increased. This can be explained by not¬ 
ing that, as the node particles align to the external field, 
the chains attached to the node are effectively shortened. 
The entropic effect of this shortening is more significant 
for a short chain than for a long one. While a long chain 
can slightly uncoil to accommodate for the motion of its 
ends, a short chain will be forced into a rather straight 
configuration, which is entropically very unfavorable. 

The influence of the network geometry (diamond cubic 
vs. simple cubic) on the amount of shrinkage is deter¬ 
mined by three factors. First, more chains attached to a 
node lead to a higher stress, when a field is applied, as 
more chains get stretched by the node particles’ rotation. 
Second, more chains lead to a higher density of parti- 


















cles in the system, and therefore also a higher stiffness. 
Finally, a contraction in one direction in most materi¬ 
als leads to an expansion in the perpendicular directions. 
The strength of this expansion in the direction perpen¬ 
dicular to the stress is measured by the Poisson ratio 
(Eqn. 14), which also depends on the network geometry. 

As can be seen from the plot of the gel’s equilibrium 
shape versus the field (right hand side of Fig.[^, the 
amount of shrinkage in the direction parallel to the field 
is larger for the diamond cubic geometry than for the 
simple cubic one. In the direction perpendicular to the 
field, gels based on the diamond cubic geometry expand 
slightly, while they contract when based on the simple 
cubic geometry. The relative shrinkage, i.e., the loss of 
volume due to an applied field, divided by the volume in 
the field free case is given by 


can be written down in matrix form as follows 


C = 


^yx 

^zx 


'^xy 

^yy 


^zy 


'^yz 

Czz 


(9) 


where the first index denotes the direction of the stress 
response and the second index denotes the direction of a 
strain. Then, for a given strain in x, ^ and 2 ;-direction 


e = e 


yy 

^zz 


the resulting stress is 


cr = Ce. 


( 10 ) 


( 11 ) 


SV{a) ^par(<^)^perp(<^)^ 

g ’ ® 

where /par(<^) and /perp(<a) denote the length of the gel 
parallel and perpendicular to an external magnetic field 
with a Langevin parameter of a, and Iq denotes the length 
of the gel in the field free case, which is equal in all 
Cartesian directions. The shrinkage for gels based on the 
diamond cubic and simple cubic geometries and chain 
lengths of 60 and 80 beads is shown on the left hand side 
of Fig.|^ It can be seen that the shrinkage is of similar 
magnitude for all systems considered, even though the 
actual shape of the gel (Fig.|^ depends strongly on the 
network geometry. The similarity of the shrinkage stems 
from the fact that, while the simple cubic gels shrink less 
in the direction parallel to the field than the diamond cu¬ 
bic ones, they shrink more in the perpendicular direction. 


V. ELASTIC CONSTANTS AND POISSON 
RATIO 


Hence, to describe the elastic properties of the gel in the 
linear response regime, the elasticity matrix C has to be 
determined. 

Depending on the symmetry of the system, some of 
the entries of the elastic matrix C are identical. For an 
isotropic material, there are only two free parameters, 
namely the diagonal and the off-diagonal elements of C. 
The diagonal elements describe the stress in the direction 
parallel to the strain, whereas the off-diagonal elements 
describe the strain in the perpendicular direction. We 
have 

fa b b\ 

C = \ b a b ] . (12) 

\b b a j 

This situation applies to the gel, when no external field 
is applied. 

When, on the other hand, a field is applied, a distinc¬ 
tion needs to be made between the direction parallel to 
the field and the two Cartesian directions perpendicular 
to it. When the external magnetic field is parallel to the 
x-direction, the resulting elastic matrix has the from 


To gain an understanding of the deformation of mag¬ 
netic gels in external fields, it is not only necessary to 
study the underlying mechanism. Additionally, the elas¬ 
tic properties of the material need to be considered, as 
they determine, to what degree the material will deform 
under a given stress. Additionally, in an external mag¬ 
netic field, not only the shape of a magnetic gel can 
change, but also its elastic properties. Hence, in this 
section, the change in these elastic constants will be ex¬ 
amined for our two gel models. The resulting stress can 
be assumed to be linear in the strains, if the strains are 
sufficiently small. In this case, the elastic constants de¬ 
scribe what kind of stress is observed in a system, when 
a certain strain or shear is applied. When only the linear 
strain but no shear is considered, the elastic constants 


f a b b\ 

C'fieid = j c d e j . (13) 

\c e d) 

As it can be seen, there are five independent elastic con¬ 
stants. The constant a describes the stress parallel to the 
field for a strain parallel to the field, b describes the stress 
parallel to the field for a strain perpendicular to the field, 
and c describes the stress perpendicular to the field for 
a strain parallel to the field. The symbol d describes the 
stress parallel to a strain which is applied perpendicular 
to the field. Finally, e describes the stress in a direc¬ 
tion perpendicular to the field, when the strain is applied 
along the second Cartesian direction perpendicular to the 
field. 
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FIG. 6. Comparison of the elastic constants a through e as 
defined by Eqn.[^for gels constructed in the diamond cubic 
and simple cubic geometries for the field-free case {a = 0). 
It can be seen, that the absolute value of the constants for 
the simple cubic system are always larger than those for the 
diamond cubic one. This indicates a more rigid gel. 


In this paper, the elastic constants are obtained by 
numerically calculating the derivatives of the stress ten¬ 
sor by running simulations at different strains close to 
the equilibrium strain, measuring the stress, and fitting 
a linear model to the resulting data. The technical de¬ 
tails of this procedure can be found in the supplementary 
information. 


The elastic constants are obtained for gels based on 
both, the simple cubic and the diamond cubic network 
geometries with a chain length of 60 beads. Measure¬ 
ments are performed for the field-free case as well as for 
an external field of a = 20. Even when no field is applied, 
no additional constraints are applied to the elasticity ma¬ 
trix. 


In Fig.[^ the elastic constants are compared for gels in 
the simple cubic and diamond cubic geometries. It can 
be seen that gels in the simple cubic geometry are signif¬ 
icantly stiffer than those in the diamond cubic one. In 
particular, the elastic constants which describe a stress 
occurring parallel to a strain (a and d in Eqn. 13) are ap¬ 
proximately four times larger. The off-diagonal elements 
of the elasticity matrix (6, c, and e) are only higher by 
10 to 25% in the simple cubic geometry. 


In Fig.[^ the elasticity of a gel in an external field 
{a = 20) is compared to that of a gel with no applied 
field. Data is presented for the DC and SC geometries, 
respectively. We observe that all ratios are larger than 
unity. Hence, the material is more stiff, when an exter- 


Diamond cubic 



Simple cubic 

FIG. 7. Gomparison of the elastic constants a through e 
(Eqn. p^ of our model gel for the cases of no external field 
and an external field of o = 20. Results for the diamond 
cubic geometry are shown at the top, while results for the 
simple cubic geometry are shown at the bottom. 


nal magnetic field is applied. This can be understood by 
noting that the gel shrinks in a field, and therefore the 
density of particles in the system is higher. It is notable 
that in the simple cubic geometry, the off-diagonal ele¬ 
ments (6, c and e) are increased more strongly than the 
diagonal elements (a and d), when an external field is 
applied. 

From the measurement of the elasticity matrix, it is 
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also possible to calculate the Poisson ratio for the inves¬ 
tigated systems. For a prescribed strain in one direction, 
these ratios measure the resulting strain in the perpen¬ 
dicular directions. In the case of a system with one spa- 
cial direction, namely the field direction, there are three 
different Poisson ratios: a prescribed strain in the field 
direction and a resulting strain in a non-field direction, a 
prescribed strain in a non-field direction and a resulting 
strain in the field direction, and a prescribed strain in a 
non-field direction and a resulting strain in the second 
non-field direction. The subscript / and n will be used 
to denote field and non-field directions, respectively. The 
first letter denotes the direction of the prescribed strain, 
whereas the second letter denotes the direction of the re¬ 
sponse. The Poisson ratio is positive, when a prescribed 
expansion leads to a contraction in the response direction. 
The field is assumed to be applied in x-direction. To ob¬ 
tain the Poisson ratio pfn for a prescribed strain in field 
direction and a resulting strain in a non-field direction, 
we apply the strain 


F— I Pfn I (14) 

V -Pfn ) 

to the system. Note that for a prescribed strain of unity, 
the resulting contraction in the perpendicular directions 
is the Poisson ratio. Now, to find the Poisson ratio, we 
require that for the given e, the stress in the response 
directions has to be zero: 



where C is the elasticity matrix. The stress occurring 
in the direction of the prescribed strain is irrelevant, so 
X can take an arbitrary value. In other words, only the 
second and third row of the equation have to be solved. 
Using Eqn.[^for the elasticity matrix C, we find 

The remaining two Poisson ratios, Pnf and Pnm can be 
obtained similarly. They are 

a^b a + e 

Based on these equations, the Poisson ratios for the 
gels can be calculated using the elastic constants shown 
in Fig.The resulting ratios are shown in FigValues 
are shown for networks based on the diamond cubic and 
simple cubic geometry, and for both, the field-free case 
and an external field oi a = 20. In the field-free case. 



FIG. 8. Poisson ratios p/n, Pnf and pnn for gel samples based 
on the diamond cubic (DC) and simple cubic (SC) geometries 
for external magnetic fields of o = 0 and a = 20. In the field- 
free case, all Poisson ratios for a given geometry are expected 
to be the same. The visible deviations are due to statistical 
measurement errors. It can be see that the direction of the 
prescribed strain as well as the direction of the measured re¬ 
sponse do not significantly influence the Poisson ratio. The 
Poisson ratio does, however, strongly depend on the network 
geometry. In the diamond cubic case, a much higher value is 
observed. 


the system is isotropic, and thus the three Poisson ratios 
Pfni Pnf 9.nd Pnn cxpccted to bc equal. In the figure, 
however, a deviation of about 20% is visible due to sta¬ 
tistical errors. It can be seen that the choice of directions 
for the prescribed strain and the observed response (field 
or non-field direction) does not have a strong influence 
on the Poisson ratio. The choice of network geometry - 
diamond cubic or simple cubic -, on the other hand, has 
a strong influence on the Poisson ratio: a much higher 
value is observed for the diamond cubic geometry. The 
Poisson ratios also increase, when a strong magnetic field 
is applied. 

In summary, a procedure has been described to obtain 
the elastic constants by fitting a linear elastic model to a 
set of stress-strain measurements. The elastic constants 
are found to be influenced mostly by the choice of net¬ 
work topology. In the diamond cubic case, the gel is much 
softer but has much larger Poisson ratios compared to the 
simple cubic case. The elastic constants turn out not to 
be strongly anisotropic, even when an external field is ap¬ 
plied. The differences in Poisson ratio help to explain the 
differences in the deformation behaviour found for differ¬ 
ent network geometries (Fig.[^. In particular, the low 
shrinkage in the direction perpendicular to the field for 
the gel based on the diamond cubic geometry is related to 
the high Poisson ratio for this system: Due to the shrink¬ 
age in the direction parallel to the field, an expansion 
occurs in the perpendicular direction. Due to the high 
value of the Poisson ratio for this geometry, this compen¬ 
sates for the shrinkage which would otherwise occur due 
to the rolling-up of polymer chains around the magnetic 
particles. 
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From the plot, it can be seen that the magnetization of 
the gels is always lower than the corresponding Langevin 
magnetization. This is due to the fact that the polymer 
chains get strained, when the magnetic particles rotate 
to align their dipole moments to the external field. The 
strain on the polymer chains causes a stress on them, 
which in turn creates a torque on the magnetic particles, 
which counteracts the magnetic field. This leads to a re¬ 
duction in magnetization. The loss of entropy per chain, 
due to its stretching, increases for shorter chains. Thus, a 
larger stress occurs and the magnetization is lower for the 
sample with shorter chains. Also, the torque on the mag¬ 
netic particles is larger, when more chains are attached 
to it. Hence, the magnetization is also lower for the gel 
based on the simple cubic geometry, in which six chains 
are connected to each node particle, compared to four 
chains attached in the diamond cubic geometry. 


FIG. 9. Magnetization curve for gels of the diamond cu¬ 
bic and simple cubic geometries for chain lengths of 60 and 
80 beads, respectively. The solid line shows the magnetization 
curve expected for non-interacting dipoles in three dimensions 
(Langevin curve). It can be seen that the magnetization for 
all gels considered is below that of non-interacting dipoles. As 
soon as the field aligns the magnetic node particles of the poly¬ 
mer network, a stress is created on the polymer chains. Due 
to this stress, the polymer chains exert a torque on the mag¬ 
netic particle which acts against the alignment to the external 
field. The torque counteracting the magnetization is larger, 
if the chains are shorter, and if there are more chains. Thus, 
a lower magnetization is observed for the case of 60 beads as 
well as the case of the simple cubic network geometry. 


VI. MAGNETIC RESPONSE 


In the model considered here, through the binding of 
the polymer chains to specific spots on the surface of 
the magnetic nanoparticles cross-linking the network, a 
coupling is created between the magnetic particles’ ori¬ 
entation and the polymer matrix. This coupling is the 
basis of the deformation mechanism for this kind of mag¬ 
netic gels. By means of the gel’s magnetic response to an 
external magnetic field, further insights into the details 
of the deformation mechanism can be gained. In Fig.[^ 
the magnetization curve M{a) is shown for gels based 
on the diamond cubic and simple cubic geometries with 
chain lengths of 60 and 80 beads. The curve describes 
the component parallel to the magnetic field of the sum 
of all magnetic moments in the system for a given field 
expressed in terms of the Langevin parameter (Eqn.[^. 
Also shown is the Langevin law, which is the result ex¬ 
pected for non-interacting dipoles, given by 


Mi{a) 


cosh a 1 
sinh a a ’ 


(18) 


VII. SUMMARY 


In this paper, we presented a detailed study of the de¬ 
formation mechanism in a magnetic gel which is cross- 
linked by magnetic nanoparticles. The deformation 
mechanism described arises from a direct coupling of the 
orientational degree of freedom of the magnetic moments 
to the polymer chains. When an external field is applied, 
the magnetic nanoparticles rotate and exert a stress on 
the polymer chains attached to them. This, in turn leads 
to a contraction of the matrix. The model, which we 
study by means of coarse-grained molecular dynamics, is 
inspired by an experimental system discussed in Ref. [? ]. 
In the three dimensional case, which is the focus of this 
report, an anisotropic deformation is observed. While 
there is a strong shrinkage in the direction parallel to the 
field, the shrinkage in the perpendicular directions is ei¬ 
ther small or not present at all depending on the network 
topology. 

A comparison was made between two network topolo¬ 
gies, one with four chains and another one with six chains 
connected to a node. In these topologies, the nodes are 
arranged in a diamond cubic and simple cubic geometry, 
respectively. An increase of the number of chains con¬ 
nected to a node has two effects. On the one hand, it 
increases the stress on the gel caused by an alignment 
of the magnetic particles to the external field. This, on 
its own, should lead to an increased shrinkage of the gel. 
However, at the same time the higher number of chains 
results in a higher particle density in the model gel, which 
should lead to a lower shrinkage. For the systems studied 
here, there is a stronger shrinkage, when four chains are 
attached to a node. As this result is presumably highly 
dependent of the details of the model, it is conceivable 
that the trade-off between the two mentioned trends can 
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lead to a different behaviour for a system with, for in¬ 
stance, significantly lower overall density. Elastic prop¬ 
erties were studied by fitting a linear elasticity model to 
a set of strain measurement. We found a significantly 
higher stiffness for the more strongly cross-linked gels in 
the simple cubic structure, as well as an increase of the 
stiffness when a field is applied. 

The magnetic response of the model gel which is cross- 
linked by magnetic particles is strongly influenced by 
the coupling between the orientation of the magnetic 
nanoparticles and the polymers. As the alignment of the 
magnetic particles to an external field exerts a stress on 
the polymers, there is an additional energy penalty for 
this alignment. In consequence, the magnetic response of 
the gel is below that of non-interacting magnetic parti¬ 
cles. 

In summary, a detailed study of the deformation, elas¬ 
ticity and magnetic response of a particle-cross-linked fer- 
rogel was performed using simulations. Based on the re¬ 
sults, we have shown that the overall deformational re¬ 
sponse is determined by an interplay between the gel’s 
degree of cross-linking and its elasticity. In the future, we 
plan to extend the model by introducing randomness into 
the chain lengths and the connectivities. Additionally, in 
experimental systems, different deformation mechanisms 
might occur at the same time. This would be the case for 
gels, which are cross-linked by magnetic particles with a 
very high magnetic moment and at a high density. 
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SUPPLEMENTARY INFORMATION 


Technical details of the determination of the elastic 
constants 


In the paper, the elastic constants are obtained by nu¬ 
merically calculating the derivatives of the stress tensor 
by running simulations at different strains close to the 
equilibrium strain, measuring the stress, and fitting a lin¬ 
ear model to the resulting data. 

To obtain the data, the stress is measured in two series, 
namely 

• a set of simulations in which one strain value is ap¬ 
plied to the Cartesian axis parallel to the field, and 
another strain is applied to both directions perpen¬ 
dicular to the field 


• a set of simulations, in which the Cartesian axis 
parallel to the field and one of the axes perpendic¬ 
ular to the field are unstrained, but the other axis 
perpendicular to the field is strained. 

The relative strains applied are ±3% and ±5% as well as 
0, resulting in 30 strain configurations. For each strain, 
eight simulations were performed for averaging. During 
the simulations, stress measurements are taken every four 
time steps. There are approximately 200 000 samples per 
simulation in the simple cubic geometry and 400 000 
samples for the diamond cubic geometry. The remain¬ 
ing simulation parameters have been described in Sec.[TT| 
The size of the unstrained system for these simulations is 
taken from the iterative procedure described in Sec.[TV| 
The equilibrium size (Fig.[^ is characterized by two pa¬ 
rameters, /par and /per, which denote the length of the 
simulation box in the directions parallel and perpendicu¬ 
lar to the applied magnetic field. A length of 1 refers to 
a simulation box, in which the polymer chains between 
the nodes of the network are in a straight line configu¬ 
ration and the spacing between neighboring beads equals 
the Lennard-Jones a. 

The error bar on these measurements of the equilibrium 
size for a given field is significant. However, these values 
need to be known to calculate relative strains in terms 
of the equilibrium size. To circumvent this problem, the 
equilibrium size is allowed to vary during the fit. This 
way, an independent estimate for the equilibrium swelling 
can be determined from the measurements of the stress- 
strain relation. 

In total, there are seven fit parameters 


P — (^par, ^per, C, d, e). 


(19) 


where /par and /per are the equilibrium size of the system, 
as described above, and a through e are the five inde¬ 
pendent entries of the elastic matrix as in Eqn. 13 The 
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target function of the fit is given by 

= ( 20 ) 

k 

where the sum is over all N strains for which a stress 
has been measured, denotes the stress measured in 
the simulation and denotes the stress calculated using 
the current estimate of the fit parameters p. That one 
is obtained by constructing the elasticity matrix C from 
the current estimates of the fit parameters a through e. 
Then the relative strain Ck for a given simulation k is 
calculated from the ratio of its actual simulation box and 
the current estimates of the equilibrium box size /par and 
^perp- The calculated stress for this simulation is then 

^k{p) =C{p)ekip). ( 21 ) 

To perform the actual fit, the target function, Eqn.[20| 
is minimized using the Lbfgs_b constraint minimization 
of the scientific python library. The entries of the elastic 
matrix are allowed to assume a value between —1 and 1. 
The equilibrium size is allowed to deviate by 0.02 from the 
estimate obtained from the iterative procedure (Fig.|^. 
In Fig.[^ the measured and calculated stresses are shown 
for the simulations in which the strain on the two axes 
perpendicular to the field is equal. It can be seen that 
there is qualitative agreement between measured and cal¬ 
culated stresses, but that there are some outliers. There 
are two causes for this, namely a large statistical error in 
the measured stress, and fitting errors. 

The quality of the fit can then be evaluated by verify¬ 
ing that those components of the elasticity matrix which 
should have identical values due to the system symmetry 
actually are similar. 

To asses the quality of the fits, the following checks are 
preformed 

• The equilibrium sizes as determined by the fit agree 
within approximately two percent with those ob¬ 
tained from the iterative procedure. 


• The entries of the elastic matrix are negative. When 
applying a positive strain, i.e., when expanding the 
system, the stress is negative. 

• When there is no magnetic field applied 

— The equilibrium size is the same in along all 
Cartesian axes (/par ~ /perp)- 

— Because of the isotropy in the field free case, 
only two values in the elastic matrix can be 
chosen freely, namely the value of the diagonal 
elements and the value of the off-diagonal ele¬ 
ments. Even when we fit the full linear model, 
_ these constants need to assume similar values, 


geo. 
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(a) 

ib) 

(c) 

id) 

(e) 

DC 

60 

0 

0.277 

0.276 

-0.282 

-0.123 

-0.109 

-0.310 

-0.101 

DC 

60 

20 

0.215 

0.276 

-0.408 

-0.205 

-0.190 

-0.522 

-0.230 

SC 

60 

0 

0.313 

0.311 

-1.186 

-0.134 

-0.140 

-1.196 

-0.131 

sc 

60 

20 

0.269 

0.301 

-1.623 

-0.300 

-0.308 

-1.763 

-0.287 


TABLE I. Elastic constants for a gel sample in the diamond 
cubic (DC) and simple cubic (SC) network geometries for a 
chain length of Ic = 60 and external magnetic fields 0 = 0 and 
o = 20. /par and /per denote the equilibrium size of the gel in 
the directions parallel and perpendicular to the applied field. 
The elastic constants Cij denote the stress observed in the i- 
direction when the system is strained in the j-direction. Due 
to the fact that the two directions perpendicular to the field 
are indistinguishable, some of the constants are equal and are 
shown only once. The letters in brackets refer to the symbols 
in Eqn.p^ 

if the fitting result is to be physically meaning¬ 
ful. I.e., in Eqn,p!3l a ^ d as well as 6 ^ c « e. 
This is the case up to an accuracy of approxi¬ 
mately ten percent for the diamond cubic ge¬ 
ometry and approximately three percent for 
the simple cubic geometry. 

The resulting elastic constants are summarized in table 
[T[ they are also shown in plots through in the paper. 
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FIG. 10. Plot of the measured stress (black arrows) and stresses calculated based on the htted elasticity matrix (red arrows). 
Results are shown for conhgurations which have a strain applied to the axis parallel to the magnetic held and an other strain 
applied to both axes perpendicular to the held. 






